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ABSTRACT 

Forming stars emit a substantial amount of radiation into their natal environment. We use ORION, 
an adaptive mesh refinement (AMR) three-dimensional gravito-radiation-hydrodyanics code, to sim- 
ulate low-mass star formation in a turbulent molecular cloud. We compare the distribution of stellar 
masses, accretion rates, and temperatures in the cases with and without radiative transfer, and we 
demonstrate that radiative feedback has a profound effect on accretion, multiplicity, and mass by 
reducing the number of stars formed and the total rate at which gas turns into stars. We also show, 
that once star formation reaches a steady state, protostellar radiation is by far the dominant source 
of energy in the simulation, exceeding viscous dissipation and compressional heating by at least an 
order of magnitude. Calculations that omit radiative feedback from protstars significantly underesti- 
mate the gas temperature and the strength of this effect. Although heating from protostars is mainly 
confined to the protostellar cores, we find that it is sufficient to suppress disk fragmentation that 
would otherwise result in very low-mass companions or brown dwarfs. We demonstrate that the mean 
protostellar accretion rate increases with the final stellar mass so that the star formation time is only 
a weak function of mass. 

Subject headings: ISM: clouds - kinematics and dynamics- starsrformation - methods: numerical — 
hydrodynamics - turbulence - radiative transfer 



1. INTRODUCTION 

On large scales molecular clouds are generally observed 
to have limited temperature variations, a characteristic 
that results from the high efficiency of radiative cool- 
ing at typical cloud densities. Consequently, simulations 
of molecular clouds frequently assume constant gas tem- 
perature, a convenient approximation for investigations 
of gas-dynamics, turbulence, and gravitational collapse 
(Gammie et al. 2003; Bonnell et al. 2003; Li et al. 2004; 
Tilley & Pudritz 2004; Vazquez-Semadeni et al. 2008). 
However, an isothermal assumption necessarily neglects 
the influence of heating due to gas compression, accre- 
tion, and stellar sources. 

The importance of the local gas temperature to the star 
formation process is motivated analytically when con- 
sidered in combination with gravity. The characteristic 
fragmentation scale for self-graviting gas of density, p, 
and sound speed, Cg, is given by the Jeans length: 



(1) 



Thus, for lower temperatures, gas is prone to gravita- 
tional instability at lower densities. In rotating self- 
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gravitating disks, the criterion may be phrased in terms 
of the local column density E: 



Q = 



(2) 




where the onset of gravitational instability occurs as the 
Toomre parameter, Q, approaches one and Ke is the 
epicyclic frequency. Cold protostellar disks more readily 
develop spiral structure and become Toomre unstable, 
influencing protostellar accretion and driving fragmen- 
tation (Kratter et al. 2008). Overproducing low-mass 
objects or brown dwarfs in the stellar initial mass func- 
tion is one side-effect of increased fragmentation (Bate 
2009b; Krumholz et al. 2007a). 

Gas eventually becomes optically thick at densities or- 
ders of magnitude above the molecular cloud mean, and 
radiative cooling is no longer efficient. To investigate 
this transition, Masunaga et al. (1998) modeled a spher- 
ically symmetric core collapse including angle-dependent 
multi- frequency radiative transfer, resolving scales down 
to the accretion shock. They halted the calculation at 
the end of the first collapse phase, prior to the dissocia- 
tion of H2 and before protostellar feedback commences. 
The authors reported a characteristic transition density 
of ~ 10^^'^ g cm^"^ for initially 10 K gas, above which 
the temperature increased with increasing density. In 
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many turbulent simulations rudimentary heating due to 
gas compression is frequently represented using an equa- 
tion of state (Li et al. 2003; Bate et al. 2003; Bate & 
Bonnell 2005; Jappsen & Klessen 2005; Bonnell et al. 
2006; Offner et al. 2008; Clark et al. 2008; Bate 2009a). 
Although such an equation typically fits a more exact ra- 
diative transfer solution like that reported by Masunaga 
et al. (1998), it neglects the instantaneous mean free 
path, multi-dimensional effects, dust chemistry, and time 
dependence of stellar sources. In fact, the subsequent 
paper, Masunaga & Inutsuka (2000), demonstrated that 
gas temperatures may bc;eomc significantly warmer as a 
result of protostellar feedback and that the temperature 
distribution is quite sensitive to the accretion luminosity. 

To compromise between physics and computational ex- 
pense, a few hybrid methods include heating by solving 
explicit diffusion approximations, estimating the instan- 
taneous radiative cooling, or extrapolating from previ- 
ously tabulated temperatures (Stamatellos et al. 2007: 
Banerjee & Pudritz 2007; Bonnell & Rice 2008). Such 
methods are computationally cheaper and reproduce ra- 
diative heating for simple geometries. However, the suit- 
ability of these approximations is unclear for radiative 
problems involving clustered star formation in a turbu- 
lent medium, where the problem is highly non-linear, in- 
volves complex geometry, and the column density may 
not be a good indicator of the cooling rate. In addi- 
tion, many of these approaches also neglect heating by 
stellar sources, which are crucial as we will show in this 
paper. The unknown accuracy of radiation approxima- 
tions and the deficiencies in handling temporal and spa- 
tial variations motivates our use of a full radiative trans- 
fer method, albeit one based upon the gray flux-limited 
difi^usion approximation. 

Relatively few authors have pursued 3D calculations 
including radiative transfer. These authors always adopt 
the flux-limited diffusion approximation and assume that 
the radiation field is frequency-independent, i.e., gray. 
By modeling star formation with gray flux-limited dif- 
fusion (GFLD), it has been shown that a barotropic 
or pol}ntropic equation of state (EOS) can underesti- 
mate the true heating at high densities even for sim- 
ple, non-turbulent collapse problems (Boss et al. 2000; 
Whitehouse & Bate 2006). The issue of radiative feed- 
back is particularly acute for high-mass stars, which 
emit prodigious luminosities while forming (Krumholz 
2006; Krumholz et al. 2007a, 2009). To explore this 
point, Krumholz et al. (2007a) contrasted simulations of 
collapsing, turbulent high-mass cores using an isother- 
mal EOS to those using GFLD radiation transfer. The 
authors demonstrated that simulations with radiative 
transfer are able to produce a massive star formed from 
gas accretion, while barotropic or isothermal calculations 
may only produce a massive star via mergers of many 
smaller bodies. Comparisons of the temperature distri- 
bution, assuming a barotropic EOS in lieu of radiative 
transfer, showed significant underestimation of the vol- 
ume of heated gas and a much lower local maximum gas 
temperature. 

In the regime of low-mass star formation. Bate (2009b) 
modeled several small clusters forming low-mass stars 
with the SPH radiative transfer method developed by 
Whitehouse & Bate (2006). The author compared these 
with previous published simulations using identical ini- 
tial conditions and a barotropic EOS. The calculations 



including radiation transfer showed a substantial de- 
crease in the number of brown dwarfs from 50% of the 
number of objects to < 10%. This agrees with the pre- 
diction of Matzner & Levin (2005), who assert incorrect 
disk fragmentation may produce brown dwarfs if irra- 
diation is not included. Accretion luminosity, which is 
emitted at the protostellar surface, generates a signifi- 
cant portion of the luminosity during the early stages 
of protostar formation. Indeed, Bate (2009b) found in- 
creased heating and fewer brown dwarfs at higher resolu- 
tion but reported little difference in the final stellar dis- 
tribution for resolutions of 0.5 AU versus 5.0 AU. Since 
0.5 AU is much larger than protostellar radii, significant 
accretion heating was neglected. As Bate (2009b) also 
neglected deuterium burning, the calculations represent 
a lower limit on the effects of radiative feedback. 

In this paper, we model the formation of low-mass stars 
in a turbulent molecular cloud including GFLD using the 
ORION adaptive mesh refinement (AMR) code. We ad- 
dress the issue of radiative feedback, including all the im- 
portant energy sources, and how it influences low-mass 
star formation. Our study differs from previous work in 
that we use source terms to account for accretion lumi- 
nosity down to the stellar surface and include a stellar 
evolutionary model (Tan & McKee 2004) . We contrast a 
GFLD simulation to one without radiative transfer. We 
also perform a less time evolved calculation with resolu- 
tion eight times higher to characterize the dependence of 
the solutions with ajid without radiative transfer on reso- 
lution. We describe our method in §2. In §3, we compare 
and contrast the four simulations. We discuss caveats to 
our method in §4 and summarize our conclusions in §5. 
Comparisons to observations will appear in a subsequent 
paper. 

2. METHODOLOGY AND INITIAL CONDITIONS 

2.1. Numerical Methods 

For the purpose of comparison, we perform two calcula- 
tions with identical resolutions and characteristic param- 
eters. The first, which we denote RT, includes radiative 
transfer and feedback from stellar sources. The second, 
henceforth NRT, uses an EOS to describe the thermal 
evolution of the gas. We perform both simulations us- 
ing the parallel AMR code, ORION. ORION utihzes a 
conservative second order Godunov scheme to solve the 
equations of compressible gas dynamics (Truelove et al. 
1998; Klein 1999): 



|+V.(pv)=0, 



dt 



+ V • (pvv) = -VP - pV(l), 



(3) 
(4) 



djpe) 
dt 



V • [{pe + P)v] = pvVcj) - KRp{4TrB - cS),(5) 



where p, P, and v are the fluid density, pressure, and 
velocity, respectively. The total fluid energy is given by 
e = l/2/c)v^ -h P/(7 — 1), where 7 = 5/3 for a monatomic 
ideal gas^ . The total radiation energy density is denoted 
by and B is the Planck emission function. ORION 
solves the Poisson equation for the gravitational poten- 
tial 4> using multi-level elliptic solvers with multi-grid 
iteration: 

V20 = 47rG[p + ^m„5(x-x„)], (6) 
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where m„ and x„ arc the mass and position of the nth 
star. 

ORION solves the non-equihbrium flux-hmited dif- 
fusion equation using a parabolic solver with multi- 
grid iteration to determine the radiation energy density 
(Krumholz et al. 2007b): 

^-V-{—VE) = «pp(47rS-c£;)+^L„W^(x-x„), 
at 



(7) 

where kr and Kp are the Rosseland and Planck dust 
opacities, and L„ is the luminosity of the nth star. W 
is a weighting function that determines the addition of 
the stellar luminosity to the grid (see Appendix A for 
details of the star particle algorithm). The flux-limitcr is 
given by A = -^(cothi? — -^), where R = \VE/{KjipE)\ 
(Levermorc & Pomraning 1981). 

We assume that the dust grains and gas are ther- 
mally well-coupled, an approximation we discuss further 
in Section 4.2. We obtain the dust opacities from a lin- 
ear fit of the Pollack et al. (1994) dust model, which 
includes grains composed of silicates, trolitcs, metallic 
iron, organics, and H2O ices. For gas temperatures 
10 < Tg < 350 K, the linear fit is given by: 

(8) 
(9) 
0.49 



KK = 0.1 + 4.4(Tg/350) cm2 g-i, 
k;p = 0.3 -I- 7.0(Tg/375) cm^ g-i. 



0.23 cm^ g ^ and Kp 



These fits give 

cm^ g~^ at the minimum simulation temperature, 10 
K. Work by Semenov et al. (2003) explores the efl:ect 
of dust composition, porosity, and iron content on the 
Planck and Rosseland average opacities. For the differ- 
ent models, they find a spread of more than an order of 
magnitude in the opacity at 10 K. The simplest model, 
based upon the assumption that dust grains are homoge- 
nous spheres, produces the lowest value for the Rosseland 
opacity, kr ~ 0.02 cm^ g~^, while the most porous and 
non-homogenous grain models produce 10 K opacities as 
large as kr ~ 0.7 cm^ g^^. For temperatures above 100 
K, the different dust models are more converged and the 
opacities arc generally within a factor of 2. As a result, 
the temperature range from 10 K to 100 K is the most 
sensitive to dust assumptions. In this range, homogenous 
models increase roughly quadratically with temperature, 
while fluffier grain models increase linearly. Our opacity 
fits are then close to the mean value of kr = 0.16 for 
the Semenov et al. (2003) models, although this value is 
more representative of porous and aggregate grains. As 
a result, our dust model is reasonable for the higher den- 
sity regions of n > 10'' cm^ typical of protostellar cores, 
but we may overestimate the dust opacity in the lower 
density cold gas by as much as a factor of 10. 

In studies of low-mass star formation, it is reasonable 
to neglect pressure exerted by the radiation field on the 
dust and gas. This is because the advection of radia- 
tion enthalpy is small compared to the rate the radiation 
diffuses through the gas: 



V • ( -^ve) 



«1, 



(10) 



^ Most of the volume of the domain is too eold to excite any of 
the H2 rotational or vibrational degrees of freedom, and thus the 
gas acts as if it were monatomic. 



where R2 — \ + X^R^ is the Eddington factor. 

Without radiative transfer, the energy exchange term 
in (5) disappears, and we close the system of equations 
with a barotropic EOS for the gas pressure: 



P 



PCs 



P \^ 2 



(11) 



where Cg — (k^T j nY^"^ is the isothermal sound speed, 
7 = 5/3, the average molecular weight /x = 2.33mH, and 
the critical density, = 2x 10~^^ g cm~^. The value of /z 
reflects an assumed gas composition of rifjc = O.lnu- The 
critical density determines the transition from isothermal 
to adiabatic regimes, and we adopt a value to agree with 
the full angle-dependent ID radiation-hydrodynamic cal- 
culation by Masunaga et al. (1998) that agrees with the 
collapse solution prior to H2 dissociation. 

We use the Truelove criterion to determine the addi- 
tion of new AMR grids so that the gas density in the 
calculations always satisfies: 



P< Pi = 



G{Axif 



(12) 



where Ax; is the cell size on level I, and we adopt a Jeans 
number of J = 0.25 (Truelove et al. 1997). In the case 
with radiative transfer, it is important to adequately re- 
solve spatial gradients in the radiation field. Radiation 
gradients are primarily associated with collapsing regions 
hosting a star but are not covered by the Jeans gravita- 
tional criterion. We find that we adequately resolve the 
radiation field and avoid effects such as grid imprinting 
by refining wherever VE/E > 0.25. Although the simu- 
lation box and gas behavior is periodic, we adopt Mar- 
shak boundary conditions for the radiation field. This 
allows the radiation to escape from from the box as it 
would from a molecular cloud. 

We insert sink, or star, particles in regions of the flow 
that have exceeded the Jeans density on the maximum 
level (Krumholz et al. 2004). These particles mark col- 
lapsing regions and also represent protostellar objects. In 
the simulation with radiative transfer, the particles act 
as radiative sources, and they are endowed with a sub- 
grid stellar model. We describe the details of this model 
and its implementation in Appendix B. By construc- 
tion, stars that approach within eight cells are merged 
together. Small sink particles, such as those generated 
by disk fragmentation, tend to accrete little mass and 
frequently merge with their more substantial neighbors 
within a few orbital times. 

2.2. Initial Conditions 

We chose a characteristic 3D turbulent Mach number, 
7W=6.6, and assume that the cloud is approximately viri- 
alized: 



5a' 



GM/R 



1. 



(13) 



The initial box temperature is T = 10 K, length of 
the box L = 0.65 pc and the average density is p = 
4.46 X 10~^° g cm~^, so that the cloud satisfies the ob- 
served linewidth-size relation (Solomon et al. 1987; Heyer 
& Brunt 2004). The total box mass is 185 Mq. 

To obtain the initial turbulent conditions, we begin 
without self-gravity and apply velocity perturbations to 
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Fig. 1. — Log gas column density of the RT (left) and NRT (right) simulations at 0.0, 0.25, 0.5, 0.75, 1.0 tg. The log density weighted 
gas temperature for the RT is shown in the center. The color scale for the column density ranges from 10-1-5 - 10°-^ g cm and 10 — 50 
K for the gas temperature. Animations of the left and right panels, as well as color figures, are included in the online version. 
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an initially constant density field using the method de- 
scribed in Mac Low (1999). These perturbations corre- 
spond to a Gaussian random field with flat power spec- 
trum in the range 1 < k < 2 where k = fcphys-^/27r is 
the normalized wavenumber. At the end of three cloud 
crossing times, the turbulence follows a Burgers power 
spectrum, P{k) cx as expected for hydrodynamic 

systems of supersonic shocks. We denote this time t = 0. 
We then turn on gravity and follow the subsequent grav- 
itational collapse for one global freefall time: 

where p is the mean box density. We continue turbulent 
driving in the simulations, using a constant energy injec- 
tion rate to ensure that the turbulence does not decay 
and the cloud maintains approximate energy equiparti- 
tion. 

The calculations have a 256"^ base grid with 4 levels of 
factors of 2 in grid refinement, giving an effective resolu- 
tion of 4096^, where Ax^ = 32 AU. In §3.2, we describe 
the results of a high-resolution core study using 7 levels 
of refinement for an effective resolution of 65, 536^ and 
minimum cell size Axj = 4 AU. Generally, the calcula- 
tions run on 128-256 CPUs. 

3. RESULTS 

3.1. Radiative Transfer and Non-Radiative Transfer 
Comparison 

In order to quantify the effects of radiative feedback 
on low-mass star formation, we compare a simulation in- 
cluding radiative transfer with a non-radiative one using 
an EOS. The latter simulation is essentially isothermal 
throughout since the highest density allowed by the Tru- 
elove criteria at the fiducial maximum AMR level corre- 
sponds to p ~ 5 X 10^^^ g cm^'^. With the adopted EOS, 
gas of this density is not heated above 11 K. 

Images of the two simulations at different times are 
shown in Figure 1. Although the simulations use iden- 
tical forcing patterns applied at the same Mach num- 
ber, the details of the turbulence differ as the two cal- 
culations have slightly different time steps and turbulent 
decay rates. Both calculations begin at t = with a 
centrally condensed region that forms the first stars, an 
imprint of the large wavenumber driving. Once gravity is 
turned on, we continue driving the simulations with the 
same energy injection rate, yielding 3D Mach numbers 
of 7.0 and 8.6 at 1 tg for the NRT and RT calculations, 
respectively. Because gravitational collapse causes non- 
turbulent velocity motions, we chose to fix the energy 
injection rate rather than the total kinetic energy. Thus, 
the root-mean-squared gas velocity no longer exactly in- 
dicates the total turbulent energy, and the Mach number 
increases above the initial value. In Tables 1 and 2, we 
list the properties of the stars formed in each calculation 
at Its- 

3.1.1. Temperature Distribution 

At t = 0, the RT simulation is nearly isothermal and 
gas temperatures, are distributed between 10-11 K (Fig- 
ure 2). Evaluated at the mean box density, the gas is 
quite optically thin with an average optical depth though 
the box of tl = L X kr/9 = 0.65 pc x 4.46 x 10~^° g 
cm^"^ X 0.2 cm^ g^^ ~ 0.02. Since the box is so trans- 
parent, the gas cools very efficiently. Small temperature 



variations arise in the initial state due to the distribu- 
tion of strong shocks. For reference, gas compressed by 
a Mach 10 shock at 10 K will undergo net heating of < 
0.1 K during a time step. Qualitatively, the change is 
so small because the radiative cooling time is a factor of 
~ 10^ smaller than the time step. 

Under the influence of gravity, collapsing regions be- 
gin to become optically thick, where individual cells at 
the maximum simulation densities reach optical depths 
of r ~ 3 when T = 100 K. Figure 2 shows the evolu- 
tion of the gas temperature distribution over a freefall 
time. There arc three processes that result in heating. 
First, there is the direct contribution from the proto- 
stars, which we add as a source term in the radiation 
energy equation. Second there is heating due to viscous 
dissipation, which is given by: 

evis = -((t' • V) • V, (15) 

where a' is the viscous stress tensor, a' = r]{S — |/V • v) 
and S = \7v + (Vv)'^ (Landau & Lifshitz 1987). We as- 
sume that the dynamic viscosity 77 = p|v| A.T/7?.eg, where 
the Reynolds number, TZcg ~ 1, at the dissipation scale. 
However, turbulent dissipation occurs over a range of the 
smallest scales on the domain, where the largest amount 
of dissipation occurs on the size scale of a cell. Thus, we 
expect this formula to be uncertain to within a factor of 
two. Third, the net heating due to gas compression is 
given by: 

ecomp = --P(V • v); (16) 

the heating is negative (i.e., cooling occurs) in rarefac- 
tions. Figure 3 shows the heating contributions summed 
over the entire domain. At t = 0, the only source of 
heating is turbulent motions. The figure demonstrates 
that after star formation commences protostellar output 
rather than compression is responsible for the majority of 
the heated gas, and at Itg protostellar heating dominates 
by an order of magnitude relative to viscous dissipation 
and four orders of magnitude relative to gas compres- 
sion. Viscous dissipation dominates the heating prior to 
star formation. After star formation is underway, viscous 
dissipation occurs in the protostellar disks. In contrast, 
turbulent shocks then contribute very little to the total. 

Figure 2 shows the evolution of the gas temperature 
distribution over a freefall time. The amount of heated 
gas (T > 12 K) increases with the number of stellar 
sources from 0.06% of the volume for one protostar at 0.5 
tg to ~ 4% at 1 tg. The corresponding mass fractions 
of the heated gas are slightly higher at 0.3 % and 5%, 
respectively. As we have seen in the previous figure, most 
of this heating is directly related to the protostars, and 
it comprises a relatively small volume filling fraction. 

The temperature distribution as a function of distance 
from the sources is shown in Figure 4. As illustrated, 
heating is local and occurs within ^ 0.05 pc of the pro- 
tostar. Since the remainder of the cloud remains near 
10 K, additional turbulent fragmentation is not affected 
by pre-existing protostars. However, radiative feedback 
profoundly influences the evolution of the protostars, ac- 
cretion disks, and stellar multiplicity as we will show (see 
§3.1.2-3.1.3). Our temperature profiles are qualitatively 
similar to those of Masunaga & Inutsuka (2000), who 
model ID protostellar collapse with radiative transfer. 
During the formation of the low-mass protostar, Ma- 
sunaga & Inutsuka (2000) also find that heating above 
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Fig. 2. — Histogram of the gas temperatures weighted by volume fraction for RT at 0.0, 0.5, 0.75, and 1.0 tff. 



TABLE 1 

RT PROTOSTAR PROPERTIES AT 1 tff 



M (Mo) M (Moyr-1)^ Mf {M0yr-l)b M {M^yr-^r L (Lq) Age (Myr) d 



1.52 


4.2 


X 


10-9 


1.1 


X 


10-7 


8.7 X 


10-6 


7.2 


0.18 


0.45 


2.0 


X 


10-8 


3.9 


X 


10-8 


4.0 X 


10-6 


0.9 


0.11 


0.09 


1.4 


X 


lo-'^ 


1.3 


X 


10-7 


8.0 X 


10-7 


0.3 


0.11 


2.91 


8.1 


X 


10-7 


1.7 


X 


10-5 


2.9 X 


10-5 


177.5 


0.10 


0.35 


5.6 


X 


10-7 


2.0 


X 


10-7 


3.5 X 


10-6 


1.3 


0.10 


2.21 


6.0 


X 


10-7 


4.2 


X 


10-6 


2.4 X 


10-5 


45.2 


0.09 


1.54 


4.0 


X 


lo-'' 


7.5 


X 


10-6 


1.7 X 


10-6 


74.6 


0.09 


1.17 


9.8 


X 


10-** 


1.7 


X 


10-5 


1.4 X 


10-5 


69.4 


0.09 


0.43 


1.2 


X 


10-6 


2.8 


X 


10-6 


6.0 X 


10-6 


8.6 


0.09 


0.48 


3.2 


X 


10-6 


7.2 


X 


10-6 


6.9 X 


10-6 


19.4 


0.08 


0.65 


1.6 


X 


10-6 


9.9 


X 


10-6 


1.1 X 


10-5 


12.9 


0.08 


0.80 


5.7 


X 


10-6 


1.7 


X 


10-5 


1.5 X 


10-5 


67.6 


0.06 


0.33 


2.1 


X 


10-5 


2.2 


X 


10-5 


2.3 X 


10-5 


79.1 


0.02 


0.06 


4.7 


X 


10-6 


5.1 


X 


10-6 


7.4 X 


10-6 


3.9 


0.01 


0.01 


3.0 


X 


10-6 


1.1 


X 


10-5 


8.6 X 


10-6 


0.8 


0.003 



^Instantaneous accretion rate. 

''Final accretion rate averaged over the last ~ 2500 yrs. 
■^Mean accretion rate averaged over the protostar lifetime. 
''Age calculated from the time of particle formation. 
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Fig. 3. — The magnitude of the heating rate due to all stel- 
lar sources, viscous dissipation, and gas compression at the times 
shown in Figure 1. 

10 K is confined within 0.05 pc of the central source 
and that significant variation in temperature occurs as 
a function of density and time. Additional studies using 

GFLD (Whitchousc & Bate 2006) or approximate radia- 
tive transfer methods (Stamatellos et al. 2007; Forgan 
et al. 2009) find similar heating beyond that expected 
from a barotropic EOS. 

Due to temperature variation with both density and 
time, wc find that the gas temperature is poorly repre- 
sented by a single EOS with characteristic critical density 
and 7. Figure 5 shows the distribution of cell tempera- 
tures as a function of cell density. For reference, we also 
plot our fiducial EOS for the NRT simulation as well as 
the EOS presented by Larson (2005). We find that many 
cells at lower densities are heated due to close proximity 
with a source. In fact, for both the EOS described in §2, 
which only includes the heating due to gas compression, 
and the Larson (2005) EOS, none of the cells are pre- 
dicted to heat much above the initial 10 K temperature. 

Nonetheless, at any given time a representative EOS 
can be formulated by fitting the mean grid cell temper- 
ature binned as a function of density. Figure 5 shows a 
least-squares fit of the temperature data for two different 
times. The magnitude of the error bars is given by the 
standard deviation of the temperatures in each density 
bin. Because such an equation fits the average temper- 
ature, there is necessarily a large scatter as illustrated 
by the error bars. The two fits return different effec- 
tive critical densities and gamma values. Thus, a single 
EOS results in a large fraction of cells unavoidably at the 
wrong temperature. 

Since accretion luminosity is predominantly emitted at 
the stellar surface, a low simulation resolution, when not 
augmented for the missing source contribution, can sig- 
nificantly neglect a large part of the heating (e.g.. Bate 
2009b). Typical pre-main sequence protostellar radii are 
expected to range from 3-5 Rq for low-mass stars (Palla 
& Stabler 1993; Robitaille et al. 2006). Thus, the tem- 
perature at a distance, r, from an emitting source, L*, is 
given by: 



1/4 



(17) 



where ctb is the Stefan-Boltzman constant, and the gas 
distribution is assumed to be spherically symmetric. 



Fig. 4. — The gas temperature as a function of distance from 
the source for all sources in the RT simulation at 1.0 tff. The line 
indicates a line with T oc r~^^^. 

Then the difference in accretion luminosity for a sim- 
ulation with minimum resolution of Rj-cs = 0.5 AU ver- 
sus a simulation resolving down to the stellar surface at 
R^ = 5 Rq is given by: 



. , Gmm 
AL=——x 

JTLrpa 



1 



Gmrh 



X (20). (18) 



Thus, the actual accretion luminosity at the higher res- 
olution is a factor of 20 larger! Since we adopt a stellar 
model to calculate the protostellar radii sclf-consistently, 
wc include the entire accretion luminosity contribution 
down to the stellar surface in our simulations. From 
(18), the difference in luminosity corresponds to a fac- 
tor of (20)^/"' or ^ 2 underestimation of the gas temper- 
ature. Nonetheless, this estimate is conservative since 
it docs not include the additional luminosity emitted by 
the protostar, which may become significant during the 
Class II and late Class I phase. Thus, we expect the sim- 
ulation of Bate (2009b) may overestimate the extent of 
small scale fragmentation and BDs formed in disks. 

3.1.2. Stellar Mass Distribution 

The large temperature range in the RT simulation has 
a profound effect on the stellar mass distribution. Figure 
6 depicts the total mass of the star-disk systems in each 
simulation, were we define the surrounding disk as cells 
with p > 5 X 10~^^ g cm~^. We find that this cutoff 
selects gas within a few hundred AU of the protostars, 
visually identified with the disk, while excluding the en- 
velope gas. Increased thermal support in the protostellar 
disk acts to suppress disk instability and secondary frag- 
mentation in the core. In contrast, protostellar disks in 
the NRT calculation suffer high rates of fragmentation. 
Most of these small fragments are almost immediately ac- 
creted by the central protostar, driving temporarily large 
accretion rates onto the central source. If we define the 
star formation rate per freefall time as 



SFR// = 



M/tff' 



(19) 



(Krumholz & McKce 2005) then the total star formation 
rate in the NRT simulation is 13% versus 7% in the RT 
simulation. Thus, the RT SFRg is almost haff the NRT 
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TABLE 2 

NRT PROTOSTAR PROPERTIES AT 1 tg 



M (Mq) M (Moyr-l)'' Mf (Mgyr-l)'' M (MQyr-i)= Age (Myr) ^ 
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''Instantaneous accretion rate. 

''Final accretion rate averaged over the last ~ 2500 yrs. 
"^Mean accretion rate averaged over the protostar lifetime. 
"^Age calculated from the time of particle formation. 



value and agrees better with observations (Krumholz & 
Tan 2007). Since the simulations have the same numer- 
ical resolution, thermal physics must be directly respon- 
sible. In the RT simulation, cores containing protostars 
experience radiative feedback that slows collapse and ac- 
cretion. 

Due to the small number statistics, we do not directly 

compare with the shape of the observed initial mass func- 
tion. Accurate comparison is also problematic because 
many of the late forming protostars arc still actively ac- 
creting. As shown in Table 1, by Its in the RT case, 
about a third, or 5 of the protostars, have accretion 
rates that are at least 5 times smaller than their indi- 
vidual mean accretion rate, indicating that the main ac- 
cretion phase has ended. Adopting an efficiency factor 
of ecore = 5 to account for mass loss due to outflows 
(Matzner & McKee 2000; Alves et al. 2007; Enoch et al. 
2008), we find that the mean protostellar mass of these 
protostars is fh = 0.4 Mq, which is comparable to the 
expected mean mass of the system initial mass function 
of ~ 0.5 M0(Scalo 1986; Chabrier 2005). 

The dynamics of close bodies and embedding gas are 
difficult to accurately resolve inside a small number of 
grid cells, so we merge particles that pass within 8 
cells. Without this limit, some of the small fragments 
would dynamically interact with the central body and 
be ejected from the stellar system. These brown dwarf 
size objects are commonly prodiiccd in simulations that 
do not include a merger criterion, typically in larger num- 
bers than are observed in the stellar IMF (e.g., Bate et al. 
2003; Bate & Bonnell 2005; Bate 2009a). As a result, the 
simulation IMF only resolves wide binaries with separa- 
tions > 300 AU. 

Figure 7 shows a histogram of all created fragments in 
both simulations, including the final mass of the objects 



that are merged. Due to the low-resolution of the disks 
in the simulations, ~ 10 cells, the many small bodies 

shown in the NRT distribution indicate numerical disk 
instability rather than small bodies forming from grav- 
itational collapse. The large number of particles that 
are created in the NRT case is directly related to the 
nearly isothermal EOS. Gravitational instability in disks 
results in filamentary spiral arms. If the gas is isother- 
mal, the filaments undergo indefinite collapse irrespective 
of the numerical resolution (Truelove et al. 1997, 1998; 
Larson 2005; Inutsuka & Miyama 1992). In a numeri- 
cal calculation, this means that all the cells along a fil- 
ament exceed the Jeans criterion nearly simultaneously 
and trigger refinement. Once the maximum refinement 
level is reached, sink particles arc introduced in cells with 
densities violating the Truelove criterion (Truelove et al. 
1997). Since our sink particle algorithm is formulated 
to represent a collapsing sphere, it is not well suited to 
filament collapse. Kratter et al. (2009, in preparation) 
have addressed this issue in thcnr pr(xlominantly isother- 
mal simulations by transitioning from an isothermal to 
an adiabatic EOS once the density reaches a factor of 
four below the density at which sink particles are cre- 
ated. This has the effect of forcing filaments to fragment 
into quasi-spherical blobs prior to sink particle creation, 
thereby allowing the collapsing objects to be faithfully 
represented by point-like sink particles. At higher reso- 
lution, the barotropic nature of our EOS is invoked and 
so much of this fragmentation disappears (see Section 
3.2.2). Similarly, in radiative calculations filamentary 
collapse is halted by heating due to radiative feedback, 
so that fragmentation is described by spherical rather 
than filamentary collapse. For either representation of 
heating, although numerical fragmentation in filaments 
is restricted, physical fragmentation may yet occur. 
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Fig. 5. — The average gas temperature at 0.8 tff and 1.0 tff as a function of density. The error bars give the temperature dispersion 
in each bin. The dashed line is a least-sqares fit of Equation 11 which returns pc and 7. The dot-dashed line plots Equation 11 with the 
original parameters: Pc = lx 10~^^ g cm~^ and 7 = 1.67. The power law density-temperature relation from Larson (2005) is also plotted. 
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Fig. 6. — The distribution of masses (star + disk) for the two 
simulations at 1.0 iff- The solid and dashed-cross lines indicate 
the NRT and RT simualtions, respectively. 

The creation and fragmentation of filaments in the sim- 
ulations is a result of gravitational instability driven by 
rapid accretion. The criterion for the onset of instability 

is similar to the classic Toomrc Q < 1 condition, slightly 
modified by the non-axisymmetry of the instabilities and 
the finite scale height of the disks, which is a result of 
turbulence driven by the accretion. This sort of gravi- 
tational instability has been investigated by Kratter et 
al. (2008, 2009 in preparation), who point out that the 
presence or absence of instability depends largely on the 
accretion rate onto the disk. The rate of mass transport 
through an a disk is 



m 



''Sjdisk 



(20) 



where Q is the Toomre parameter for the disk and Cg.disk 
is the sound speed within it. Gravitational instabilities 
produce a maximum effective viscosity a ~ 1. At early 
times, wc find that the accretion rate from a core onto the 
disk forming within it can be ^ Cg core/^) where Cs,core 
is the sound speed in the core. If the sound speeds in the 
disk and core arc comparable, Cg.disk ~ Cs^corej as is the 
case in the low-resolution NRT simulation, then the disk 
can only deliver matter to the star at a rate ~ c^core/^ 



while still maintaining Q > 1. As a result matter falls 
onto the disk faster than the disk can deliver it to the 
star, and the disk mass grows, driving Q toward 1 and 
producing instability and fragmentation, as illustrated by 
the NRT simulation. Conversely, if the disk is warmed, 
either by radiation or by a switch from an isothermal 
to an adiabatic equation of state, then Cs,disk > Cs,core 
and the rate at which the disk can deliver gas to the 
star increases. If the disk is sufficiently warm then it can 
process all the incoming material while still maintaining 
Q > 1. As a result the disk docs not fragment, as is seen 
in both the low- and high-resolution RT simulations and 
in high-resolution NRT simulation. This shows that the 
fragmentation in the low-resolution NRT simulation is 
indeed numerical rather than physical in origin, and that 
it is a result of the density-dependence of the equation 
of state rather than of the resolution directly. 

This analysis also sheds light on the importance of nu- 
merical viscosity. Krumholz et al. (2004) show that in 
the inner few cells of disks, numerical viscosity can cause 
angular momentum transport at rates that correspond 
to a > 1. However, as the analysis above shows, in- 
creasing a tends to suppress fragmentation rather than 
enhance it. Wc find that fragmentation is more preva- 
lent in the low-resolution NRT simulation than the high- 
resolution one, which is exactly the opposite of what we 
would expect if numerical angular momentum transport 
were significantly influencing fragmentation. Therefore 
wc conclude that numerical angular momentum trans- 
port is not dominant in determining when fragmentation 
occurs in our simulations. 

In isothermal calculations, the issue of filamentary col- 
lapse is a problem for all sink particle methods and it 
is not unique to grid-based codes. Due to the filamen- 
tary fragmentation in the NRT case, we prefer to merge 
close particle pairs in the simulations rather than follow 
their trajectories. Note that particles are inserted with 
the mass exceeding the Truelove criterion rather than the 
net unstable mass in the violating cells. Particles created 
within a discrete bound mass typically gain size quickly. 
Most particles formed in the unstable disk regions form in 
a spiral filament and do not have significant bound mass, 
so the particle mass is tiny when they are accreted by the 
central object. However, if several small particles are ere- 



10 



atcd within the merging radius each time step around a 
particular protostar, their merging can significantly in- 
crease the instantaneous accretion rate. As illustrated by 
the figure, there are only a handful of objects that form 
and approach within a merging radius in the RT simu- 
lation, whereas the NRT simulation produces a plethora 
of such bodies. 

Bate (2009b) finds a similar reduction in protostar 
number with the addition of radiative transfer. As in our 
calculation, the final number of stars including radiative 
transfer is sufficiently small that a statistical compari- 
son with the IMF is problematic. Instead, we base our 
comparison on the mean stellar mass. Using a resolution 
of 0.5 AU, Bate (2009b) finds rh - 0.5 M©, which does 
not include outflows or any scaling factor accounting for 
their presence. Adopting a scaling factor of Ccoro = 1/3 
would produce a mean of m ^ 0.2 M0, lower than our 
RT mean mass and the mean mass of either the sys- 
tem or individual stellar initial mass function reported 
by Chabrier (2005). However, in Bate (2009b) a num- 
ber of the protostars are continuing to accrete and have 
not reached their final mass. In addition. Bate (2009b) 
demonstrates that the mean stellar mass increases as cal- 
culations approach higher resolution and include a larger 
portion of the accretion luminosity. This result is most 
likely because disk fragmentation decreases as the gas 
becomes hotter, thus increasing accretion onto primary 
objects. It is possible that if Bate (2009b) had included 
all the accretion and stellar luminosity, the mean mass 
obtained would be closer to the value we find. 

Observations suggest that BDs compose ^ 30% of the 
total population of clusters (Andersen et al. 2006). De- 
spite the merger criterion we adopt, the NRT calculation 
produces a significant number of BDs, > 30% sans scal- 
ing with ecore, resulting in a slightly lower mean mass 
than the RT run. In comparison, Bate et al. (2003) find 
that approximately half of the objects formed arc BDs. 
resulting in a mean mass of 0.1 Mq. This result per- 
sists for barotropic calculations modeling more massive 
clusters with superior resolution (Bate 2009a). Calcu- 
lations using a modified EOS that includes effects due 
to the internal energy and dissociation of H2, ionization 
state of H, and approximate dust cooling find increased 
disk fragmentation, leading to numerous BDs (Attwood 
et al. 2009). Thus, the overproduction of BDs in non- 
radiative simulations substantiates the importance of ra- 
diative transfer and feedback from protostars in accu- 
rately investigating fragmentation and the initial mass 
function. 

3.1.3. Accretion Rates 

As indicated by the Toomre criterion given by equation 
(5.2), the local gas temperature is key to the stability of 
disks, dumpiness in the disks is directly reflected in the 
variability of the protostellar accretion rate. Figure 8 
shows the accretion rates for the two first-forming proto- 
stars in each calculation as a function of time. The RT 
protostellar accretion in the left panel illustrates that 
once a protostar has accreted most of the mass in the 
core envelope, its accretion rate diminishes significantly. 
Protostars in both simulations show evidence of variable 
accretion on short timescales. However, the accretion 
bursts in the NRT simulation may vary by an order of 
magnitude, while in the RT case variability is generally 
only a few. Disk dumpiness may be magnified due to 



dynamical perturbations by nearby companions. For the 
cases shown, the RT protostar is single, while the NRT 
protostar has several companions. Similar variability to 
the NRT protostellar accretion rate is also observed by 
Schmeja & Klessen (2004). In their turbulent isothermal 
runs, Schmeja & Klessen (2004) show that the magnitude 
of the initial particle accretion rate is comparable to our 
calculations at to ^ few x 10^'"' Mq yr^^ with variability 
by factors of 5-10. However, the reported accretion rates 
appear to significantly decrease within 0.1 Myr. 

In principle, a sizable amount of the protostellar mass 
may be accreted during the periods of high accretion. 
We define a burst as an increase of 50% in the accretion 
rate over 1000 years, where mergers of another protostar 
of mass TO > 0.1 Mq are excluded. Using this metric, the 
NRT protostars accrete from 0-13 % of their mass during 
the bursts with a median of 5%. The RT protostars ac- 
crete 0.0-9 % of their mass during bursty accretion with 
a median amoimt of 1%. Thus, variable accretion is not 
significant. Our data analysis is limited by the coarse 
level time step of ~ 100 yrs, so that accretion rate vari- 
ability on shorter timescales will not be resolved in the 
analysis. For comparison, Vorobyov & Basu (2006), mod- 
eling the formation and accretion history of a protostar in 
two-dimensions, find that > 50 % of the protostellar mass 
is gained in short intense accretion bursts. In their sim- 
ulations, accretion occurs smoothly until t ~ 0.15 Myr, 
where variability on timescales < 100 yrs begins, corre- 
sponding to accretion of ^ 0.05 Mq clumps. Although 
their time resolution is finer, sampling at longer time 
increments, as in our calculation, is unlikely to miss per- 
sistent cyclical variability of four to five orders of mag- 
nitude in accretion rate. We find that when the stellar 
mass is about half the final mass, large variability in the 
RT accretion rates is rare, while it is more common in 
the NRT case. RT protostars with ages comparable to 
t ~ 0.15 Myr experience the most variable accretion oc- 
curring over 1-2 orders of magnitude. However, by this 
time, the majority of the envelope mass has been ac- 
creted and accretion rates are m ~ 10"'' M© yr~^, so 
that accreting significant mass is unlikely. 

In Vorobyov & Basu (2009), the authors demonstrate 
that simulations with a stiffer equation of state and 
warmer disk exhibit variability over at most two orders 
of magnitude. This finding is more consistent with our 
results, and it supports the differences in accretion we 
find between the NRT and RT calculations. However, 
bursty accretion due to disk instability also depends upon 
the core rotation and the rate at which mass is fed into 
the disk from the cnvelope(Vorobyov & Basu 2006: Bo- 
ley 2009). Thus, we expect that radiative effects alone 
cannot completely determine accretion behavior. Since 
the disks in our low-resolution calculations are not well 
resolved, it is possible that we may not be able to re- 
solve the disk dumpiness observed by Vorobyov & Basu 
(2006). Their innermost cell is placed at 5 AU, which is 
comparable to the cell size in our high-resolution runs, 
however, they adopt logarithm spacing to concentrate 
cells in the inner region of the disks. We note that their 
method also includes an approximate treatment of mag- 
netic fields that could influence their results and which 
we neglect in our calculations. 

Figure 9 shows that the NRT simulation exhibits 
slightly higher average accretion. Note that we sub- 
tract the accretion spikes caused by significant merg- 
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Fig. 7. — The figures show the distribution of particles formed as a function of mass for the RT(Ieft) and NRT(right) simulations (solid 
line). These include the particles that are merged, where the total particle number with final masses greater than 10~^ Mqis 23 and 251, 
respectively. The dashed lines show the distributions of stellar masses at the final time output. 
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Fig. 8. — The accretion rate, M, as a function of time for the first forming object in the RT (left) and NRT (right) simulations, 
average both simulations over 2 kyr for consistency. 



We 



ers. The mean accretion rate over the protostars hfetime 
for the final protostars is ^ 1 x 10~^ M0 yr"'^ versus 
~ 6 X 10"^ M0 yr-i for the RT run. Without the added 
thermal support from radiation feedback and with in- 
creased fragmentation, the NRT protostars accrete their 
envelope mass more quickly. However, the protostars 
in both calculations satisfy the same accretion-mass re- 
lationship, with accretion increasing approximately lin- 
early with star mass. Using a least-squares fitting tech- 
nique, we obtain power-law relationships m oc m°'^^ 
and TO oc m"-^^ for the RT and NRT data, respectively, 
which have values of 67.6 and IS.O.^We include masses 
TO > 0.1 MqIji the fit and weight the data by the ages 
of the protostars. Thus, young protostars with only a 
short accretion history are weakly weighted. As Figure 
9 shows, there is a significant amount of scatter about 
the fits. Schmeja & Klessen (2004) find a similar trend 
between the mean accretion rate and final masses for 
protostars forming in their isothermal driven turbulence 
simulations. 



The X value for the fit is given by; x 







where j/j are the age-weighted accretion rates, Xi axe the masses, 
A and B are the fit coefficients, and Cy is the standard deviation 



The apparent correlation between stellar mass and av- 
erage accretion rate occurs because protostars forming 
in more massive cores tend to be more massive and also 
have higher accretion rates. McKee & Tan (2003) derive 
a self-similar solution for the accretion rate where the 
pressure and density each have a power-law dependence 
on r, such that p oc r~^i' and P oc p^'' cx r"*^^, where 
7p = 2fcp/(2 + kp) and kp = 2/(2 - 7p). Although the 
simulated cores are not self-similar, it is possible to fit a 
power-law to the pressure of the core envelope in most 
cases. Both RT and NRT cores have exponents in the 
range A;p~0 — 5ata few thousand AU from the proto- 
star, with an average value of fcp ^ 1 or kp = 1.5. McKee 
& Tan (2003) show that the accretion rate is then: 



TO.=5.5x 10-^.^1/8^ 

Ps,corc / \ I 



•V4 



106 K cm" 



1*/ 



1 Mo, 

3(2-fcp)/[2(3-fep)] 



M0yr(^],) 



where to*/ is the final stellar mass, Ps,core is the core 
surface pressure, (jf)* and A are order unity constants de- 

of the j/i values. 
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scribing the effect of magnetic fields on accretion and the 
isothermal density profile, respectively. Since we weight 
the fit by the protostellar age, this selects for the case 
where m* ~ m*/. Assuming that Sd is roughly con- 
stant, m* oc rnj^, that is similar to the slopes produced 
by the least-squares fit. 



3.1.4. Multiplicity 

The number of stars with stellar companions is an im- 
portant observable that may directly relate to the ini- 
tial conditions of star-forming regions. Among the pop- 
ulation of field stars, most systems are single with the 
number of systems containing multiple stars increasing 
as a function of stellar mass (Lada 2006). Young pre- 
main sequence stellar populations are observed to con- 
tain more multiple systems than field stars suggesting 
that the multiplicity fraction evolves over time (Duchene 
et al. 2007). Unfortunately, the initial stellar multiplicity 
is challenging to directly measure due to the difficulty of 
resolving close pairs and limited sample sizes (Duchene 
et al. 2007). The two dominant effects influencing multi- 
plicity are fragmentation and N-body dynamics. While 
fragmentation in a collapsing core may result in multi- 
ple stars, systems with three or more bodies are dynam- 
ically unstable, causing higher-order stellar systems to 
rapidly lose members. Multiple stellar systems can also 
occur via stellar capture, a mechanism most applicable to 
high-mass stars forming in very clustered environments 
(Moeckel & Bally 2007). Goodwin & Kroupa (2005) sug- 
gest that that observed higher-order multiple systems are 
initially members of open stellar clusters rather than aris- 
ing from the fragmentation of a single core. In general, 
the number of such systems is observed to be small, with 
only 1 in every 50 systems in the field having at least 
four members (Duquennoy & Mayor 1991). 

The RT and NRT calculations present very different 
pictures of the initial stellar multiplicity. The large dif- 
ferences in temperature and fragmentation have a signifi- 
cant effect on the fractions of stars in single and multiple 
systems. As shown in Figure 10, the majority of stars 
formed in the RT calculation are single, while in the NRT 
calculation the majority of stars live in systems with 2 
or more stars. This is mainly due to continued disk frag- 
mentation rather than long-lived stable orbital systems. 
The field single star fraction (SSF), defined as the ratio 
of the number of primary stars without a stellar compan- 
ion to the total number of stellar systems, is observed to 
be - 70% (Lada 2006)^. The RT calculation produces 
an SSF of 0.8 +0.2/ - 0.4, while the NRT calculation has 
an SSF of 0.6 ±0.4, where the uncertainty is given by the 
poisson error. Due to the resolution of our calculation, 
we can only capture wide binary systems of r > 300 
AU. However, a number of protostars have undergone 
significant mergers, which we define as those in which 
the smaller mass exceeds 0.1 Mq. We find that about a 
third of the stars in the RT simulation and a tenth of stars 
in the NRT simulation have experienced significant past 
mergers. Assuming that these would have resulted in 
multiple stellar system revises the SSF values to 0.5±0.3 
and 0.6±0.4, respectively. Unfortunately, this is a very 
uncertain estimate as we have small statistics, and we 

^ The SSF does not include brown dwarfs in estimating multi- 
plicity. 
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Fig. 9. — The plot shows the distribution of average accretion 
rates (crosses) as a function of final star mass at 1 tg. The hor- 
izontal line indicates the Shu (1977) accretion rate c^/G at lOK. 
The dashed and dot-dashed lines indicate the age weighted fit of 
the average accretion rates for the RT and NRT runs, respectively. 
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Fig. 10. — The plot shows the system multiplicity for the two 
calculations, where N is the number of stellar systems, and the plot 
is normalized to the total number of systems. The multiplicity on 
the X-axis is the number of stars in each system. 

cannot know whether the systems with significant merg- 
ers would have resulted in bound or unbound systems in 
the absence of the mergers. 

3.1.5. Stellar Feedback 

Our model includes accr(rt,ion luminosity and a sub- grid 
stellar model estimating the contribution from Kelvin- 
Helmholz contraction and nuclear burning (see Appendix 
B) The stellar model includes four evolutionary stages. 
The earliest stage occurs when the protostar begins burn- 
ing deuterium within the core at a sufficient rate to main- 
tain a constant core temperature. Once the initial deu- 
terium in the core is depleted, burning occurs at the rate 
that new matter convects inwards; this is the steady core 
deuterium state. In the third stage, the star burns the 
deuterium remaining in the outer layers of the protostar. 
Finally, the star ceases contracting and reaches the zero- 
age main sequence (ZAMS). 

Figure 11 shows the luminosity as a function of time 
for three different protostars. At early times, accretion 
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Fig. 11. — The plot shows the total luminosity as a function of 
time for three stars in the RT simulation. The accretion luminosity 
contribution is shown by the dashed line, and the masses are 1.5, 
0.45, and 0.35 Mq, respectively. The bottom plot shows the total 
luminosity including all the protostars. 

dominates the luminosity, and variability in accretion is 
strongly reflected in the total luminosity. At late times, 
accretion slows and Hayashi contraction begins to make a 
substantial contribution. In general, the total luminosity 
summed over all the stars is dominated by those proto- 
stars with the highest accretion rates. For these young 
sources, the stellar luminosity is quite small in compari- 
son to the accretion luminosity. Thus, the last panel in 
Figure 11 shows that for all times, accretion luminosity 
is the main source of luminosity. 

For comparison, luminosity due to other physical pro- 
cesses such as compression and viscous dissipation is 
small (sec Figure 3). Figure 12 shows the final lumi- 
nosity as a function of source mass. The luminosity in- 
creases roughly linearly with mass but has a fairly large 
scatter. As indicated on the plot, two of the stars have 
reached the ZAMS, which was due to increased accre- 
tion resulting from significant mergers. Even in this low- 
mass stellar cluster, there are individual stars with con- 
tributions larger than the net viscous dissipation. This 
demonstrates that any heating due to viscous dissipation 
is exceeded by modest protostellar feedback. 

3.2. Resolution and Convergence 

The AMR methodology allows flexibility in both the 
depth and breadth of resolution. An insufiicient amount 
of resolution may give inaccurate results, so it is impor- 
tant to gauge the sensitivity of the result to the resolu- 
tion. The large scope of the problem and the expense 
of the radiative transfer methodology limits the depth 
or maximum resolution of our calculation, where the RT 
calculation cost is ~ 70,000 CPU hrs on 2.3 GHz quad- 
core processors. To quantify the effects of resolution on 
the solution, we run second RT and NRT calculations 
that evolve the first formed object to a resolution eight 
times higher than the overall calculation. We run these 
simulations for 0.12 ts after the formation of the proto- 
star. We adopt a fixed number of cells for the closest re- 
solved approach between two particles, so that the high- 
resolution simulations have a merging radius of 32 AU, 
a factor of eight smaller than the low-resolution cases. 

3.2.1. High-Resolution Study with Radiative Feedback 
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Fig. 12. — The plot shows the distribution of luminosities 
(crosses) in the RT simulation as a function of final star mass at 
1.0 iff. The crosses, stars, and diamonds refer to stars undergoing 
variable core deuterium burning, undergoing steady core deuterium 
burning, or reaching the zero-age main sequence. 

The high-resolution and low-resolution calculations 

both form single objects with stable, thermally supported 
disks. Figure 13 shows a comparison of the densities, 
temperatures, and radiation fields. The effective radia- 
tion temperature differs by only a few percent outside 
the inner cells of the low-resolution calculation. In both 
cases, the gas and radiation temperatures are well cou- 
pled such that Tgas — T^ad- However, the gas in the 
high resolution case is more centrally concentrated, and 
the disk radius appears smaller. At the final time, the 
high-resolution star has accreted 0.54 Mq, while the low- 
resolution case has reached 0.50 Mq. During the course 
of the run, the lower resolution case forms a few frag- 
ments in the disk, which arc almost immediately accreted 
by the primary, while in the high-resolution case, no ad- 
ditional particles arc formed. 

Figure 14 shows a comparison of the accretion and lu- 
minosity as a function of time. Accretion is generally 
smooth, and the rates are generally within a factor of 
two. The luminosity in the low-resolution run has slightly 
larger variation, but the two approach a similar value at 
later times. Although there arc deviations in the history 
between the two runs, the evolution is not significantly 
different at the higher resolution. Certainly, even higher 
resolution is preferable for investigation of disk proper- 
ties, but our main result-that radiative feedback is im- 
portant to the formation of low-mass stars-is insensitive 
to the simulation resolution. High-resolution radiation- 
hydrodynamics simulations of low-mass disks including 
irradiation confirm that such disks, with properties sim- 
ilar to ours, are stable against fragmentation (Cai et al. 
2008). Gravitational instability is expected to occur only 
in the regime where the mass of the disks is comparable 
to the stellar mass (Cai et al. 2008; Stamatellos & Whit- 
worth 2008, 2009). 

3.2.2. High- Resolution Study with a Barotropic EOS 

This higher resolution non-radiative study achieves 
maximum densities > 5 x 10~^^ g cm~^, several times 
higher than the barotropic critical density. Consequently, 
dense gas is heated to temperatures of ^ 20 — 25 K. Dur- 
ing the time we compare the non-radiative simulations, 
both the high-resolution barotropic calculation and the 
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Fig. 13. — From left to right, the images show the log density, log radiation temperature, Tr = (Br/a)^''*, and log gas temperature for 
a protostellar system at ~ 0.6 iff followed with dx = 4 AU resolution (top) and dx = 32 AU (bottom). The image is 0.03 pc on a side, 
where we denote the star position with a black cross. The colorscale ranges are given by 10~^^ — 10~^^ g cm""', 1 — 100 K, and 1 — 100 
K, respectively. 



first collapsing core in low-resolution NRT calculation 
form a similar mass primary object with protostellar disk 
(see Figure 15). However, the low-resolution NRT sys- 
tem experiences significantly more fragmentation. We 
find that the protostellar disk in the NRT case fragments 
during approximately half of the time steps, while in the 
higher resolution barotropic case fragmentation occurs 
very rarely, taking place in less than <0.1% of the time 
steps. 

Since the low-resolution NRT disks are essentially 
isothermal, we conclude that heating due to the 
barotropic approximation is largely responsible for de- 
creasing the number of fragments. In contrast, the higher 
resolution disks are heated to ^ 20 K. However, this is 
still significantly less heating than in the RT case, and 
we find that numerical instability is not suppressed com- 
pletely even with high-resolution. The radiative high- 
resolution case experiences no disk fragmentation, under- 
scoring our conclusion that radiative feedback is crucial 
to representing fragmentation or lack thereof in the star 
formation process. 

Despite different merger radii, in both non-radiative 
cases all of the fragments are eventually merged with 
the primary protostar so that the end result in both cal- 
culations is a single protostar. This suggests that the 
fragmentation taking place at low-resolution is largely 
numerical rather than physical. We emphasize that both 
significantly higher resolution than we use and additional 



physics are required to study accretion disk properties. 

3.2.3. Convergence 

The minimum breadth of resolution is determined by 
the Truelove criterion. Due to the radiation gradient re- 
finement criterion we apply to resolve the radiation field, 
at 1 tff the RT simulation has ~ 80 % more cells, gen- 
erally concentrated near the protostars, than the NRT 
calculation. This extra refinement improves the resolu- 
tion regions near protostellar sources. Inverting (5.12) 
yields an expression for the effective Jeans number for 
each cell as a function of density, resolution, and sound 
speed: 

•^cff = ,^^1/2 • (22) 

As shown in Figure 16, the RT simulation is shifted to 
lower Jeff, where the vast majority of cells in both cal- 
culations are resolved to better than Jeff = 0.1. The 
choice of base grid resolution guarantees that Jeff is typ- 
ically much smaller than J for most of the cells on the 
domain. We use a fiducial value of J = 0.25 to trigger 
additional refinement in both simulations, so no cell has 
Jeff exceeding 0.25. Cells in the highest J^s bin are ex- 
clusively found on the maximum AMR level, and they 
are generally at the highest gas densities. These cells, 
many located in the disks around the protostars, are at 
the same resolution in both calculations. Thus, the frag- 
mentation results of the RT and NRT calculations are 
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not dissimilar due to differences in effective resolution 
but are solely a result of differences in thermal physics. 

4. SIMPLIFYING ASSUMPTIONS 

These numerical calculations neglect a number of ar- 
guably crucial physical processes in low-mass star forma- 
tion. In this section, we discuss the implications for our 
results. 

4.1. Chemical Processes 
4.1.1. Dust Morphology 

Our dust model neglects the evolution of dust grains 
due to coagulation and shattering. In cold dense environ- 
ments, such as protostellar disks, the aggregation of dust 
grains may significantly increase grain sizes on timescales 
as short as 100 years (Schmitt et al. 1997; Blum et al. 
2002). Observations of Class protostars indicate sig- 
nificant evolution of the dust size distribution at average 
densities of n ~ 10^ cm~^ by the Class phase (Kwon 
et al. 2009). Since we adopt a single dust model for the 
entire domain, we are likely to either overestimate or un- 
derestimate dust grain size in different regions. 

To examine the effect of the dust model on gas tem- 
perature, we repeat the turbulent driving phase (without 
gravity) using a conservative model more typical for non- 
aggregate dust grains: 

KR = 0.015(r2/110) cm2 g-i for Tg < 110 (23) 

Kp ^ O.lO(TgVllO) cm2 g-i for Tg < 110. (24) 

Using this model, we find that shocked gas may be heated 
as high as 18 K after a crossing time. In comparison, gas 
in the fiducial case is only heated to ~11 K at the same 
densities (see Figure 2 for the temperature distribution 
due to the fiducial dust opacity model). However, the 
extent of the additional heating is quite small since only 
0.003% of the mass is heated above 11 K and thus differs 
from the fiducial case. This suggests that the simulations 
may underestimate the gas temperature in low density 
regions outside of cores {uh < 10'') where the dust dis- 
tribution is not expected to evolve due to coagulation. 
Significant discrepancy between the gas temperatures of 
the two models is mainly confined to a small number 
of cells and is mitigated by the importance of molecu- 
lar cooling in these regions, which we discuss in Section 
4.1.2. 

4.1.2. Gas Temperature at Low Density 

To simplify the dust-gas interaction, we assume that 
dust and gas are perfectly coUisionally coupled such that 
their temperatures are identical. In molecular clouds, 
there can be significant variation between the dust and 
gas temperatures. For example, dust in close proximity 
to stellar sources is radiatively heated, while in strongly 
shocked regions of the flow, dust acts as a coolant for 
compressionally heated gas. Below we will discuss both 
the regime where dust cooling dominates, i.e., Tg > T^, 
and where molecular cooling dominates, i.e., T^ > Tg. 

When the gas is shock-heated, the perfect coupling ap- 
proximation remains valid as long as the rate of energy 
transfer between the gas and dust is balanced by the 
cooling rate of the dust. The dust cooling per unit grain 
area by photon emission is: 

F(a, Td) = 4 < Q{a, T^) > aeTj, (25) 
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Fig. 14. — The plot shows the accretion rate and luminosity as a 
function of time for the first formed star in the RT calculation and 
the same object followed with da; = 4 AU resolution. Temporal 
bins of Ikyr are used. 




Fig. 15. — The images show the log density (left) and log gas tem- 
perature (right) for a NRT protostellar system at 0.5 iff followed 
with dx = 4: AU resolution (top) and dx = 32 AU (bottom). The 
image is 0.03 pc on a side, where we denote the star position with 
a black cross. The color scale ranges are given by 10-1^ - 10-1* g 
cm" and 1 — 50 K, respectively. 

where Td is the dust temperature, a is the grain size, and 
< Q(a, Td) > is the Planck-averaged emissivity (Draine 
& Lee 1984). Then for an ensemble of grains with dust 
opacity, Kp, the dust cooling is given by: 
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In equation (27) , we substitute Equation (24) for Kp and 
assume that Tg Td. The rate at which energy is trans- 
ferred from the gas to the dust is given by: 
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whore we adopt ad = 2.44 x 10^^^ cm"^ for the dust 
cross section per H nucleus (Young et aL 2004). For a 
gas temperature of 10 K the exponential term is very 
small, so we neglect it in the following equation. Equat- 
ing these expressions and solving for the gas density at 
which heating and cooling balance gives: 
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Thus, we demonstrate that the dust and gas are well 
coupled as long as the gas densities arc siifficicntly high. 

However, even in regions where the dust and gas may 
not be well coupled, molecular line cooling is important. 
For gas densities in the range rin = 10^ — 10^ cm~^, CO 
is the dominant coolant. For these densities, the cooling 
rate per H, A/nn, is approximately constant with den- 
sity at fixed temperature. To compare to the magnitude 
of dust cooling consider a 2 km s^^ shock that heats the 
gas above 100 K. The cooling rate at 100 K is given by 
n^^moi — 5 X 10~^''nH ergs cm~'^ s~^, where we adopt 
the cooling coefficient from Neufeld et al. (1995). The 
characteristic cooling time is ^ 1000 yrs at the average 
simulation density, which is approximately half the cell- 
crossing time of such a shock, implying that molecules 
cool the gas relatively efficiently. Since the shock temper- 
atures on our domain are limited by our resolution, which 
is much larger than the characteristic cooling length, 
post-shock temperatures do not surpass 20 K. In this 
regime, the dust cooling for perfect dust-gas coupling is 
at least an order of magnitude larger than the estimated 
molecular cooling. As a result, we likely under- estimate 
the temperatures in low-density strongly shocked gas in 
comparison to similar shocks in molecular clouds. 

In the regions near protostars, the perfect coupling as- 
sumption is valid provided that gas heating by dust is 
balanced by molecular cooling. This case is discussed by 
Krumholz & McKee (2008), where the authors demon- 
strate that the dust and gas temperature remain within 
a degree provided that: 
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for gas temperatures around 10 K. For higher gas tem- 
peratures around 100 K, we adopt the molecular cooling 
coefficient above and find that the dust and gas are well 
coupled provided njj exceeds ^ 2 x 10*^ cm^''. Number 
densities of this magnitude are exceeded in collapsing 
cores, so that regions near protostars are guaranteed to 
have well-coupled dust and gas. 

Thus, gas temperatures in our RT simulation are fairly 
accurate for densities larger than the average density, 
but they may be underestimated by a factor of ^ 2 in 
strong shocks when the molecular cooling rate is much 
smaller than the implemented dust cooling rate. Since 
gas heating suppresses fragmentation, our results may ac- 
tually overestimate the amount of fragmentation. Conse- 
quently, our finding that radiative feedback reduces frag- 
mentation would genca'ally be strengthened by a better 
treatment of the thermodynamics. 

4.2. Magnetic Fields 
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Fig. 16. — Histogram of the effective Jeans number, Jgff at 
1.0 tff. The solid and dashed Uncs indicate the NRT and RT 
simulations, respectively. Each histogram is normalized to the total 
number of cells. 

Observations indicate the presence of magnetic fields 
in nearby low-mass star-forming regions (Crutcher 1999). 
However, the magnitude of the fields and their impor- 
tance in the star formation process remain uncertain. 
Observations by Troland & Crutcher (2008) suggest that 
the energy contributed by magnetic fields on core scales 
is subdominant to the gravitational and turbulent ener- 
gies. On smaller scales, magnetic fields are believed to 
be associated with disk accretion and the generation of 
protostellar outflows (Shu et al. 1994; Konigl & Pudritz 
2000). 

Numerical simulations have demonstrated that the 
presence of magnetic fields may suppress disk fragmenta- 
tion by supplying additional pressure support (Machida 
et al. 2008; Price & Bate 2007, 2008). We find that the 
inclusion of radiative transfer has a similar stabilizing 
influence on disks. 

4.3. Multi-frequency Radiative Transfer 

Due to the expense of the calculations, we adopt a 
gray radiative transfer fiux-limited diffusion approxima- 
tion. By averaging over angles and frequencies to ob- 
tain the total radiation energy density, we sacrifice the 
direction and frequency information inherent in the ra- 
diation field. As discussed in Krumholz et al. (2009), 
these approximations touch on several competing effects 
that influence the radiation spectrum. Since radiation 
pressure is negligible for low-mass stars, it does not af- 
fect the gas dynamics. Instead, our main consideration 
is the extent to which radiative heating may differ for 
a more sophisticated radiative treatment. As we have 
discussed in previous sections, the gas temperature and 
corresponding thermal pressure alone have a significant 
relationship with accretion and fragmentation. 

The first effect to consider is a more exact treatment 
of dust opacity, which is strongly frequency dependent 
in the infared and increases towards lower wavelengths 
(e.g., Ossenkopf & Henning 1994). Since long- wavelength 
radiation has a lower optical depth, in a multi-frequency 
calculation the longest wavelengths would be able to es- 
cape the core. Anisotropies in the radiation field may 
also facilitate cooling. Radiative beaming, for example 
via an outflow cavity, may allow photons to escape along 
the poles (Krumholz et al. 2005). Thus, both these ef- 
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fects will likely decrease the temperature in protostellar 
cores. 

The gray radiative transfer method also assumes that 
the radiation field is thermalized, producing a Planck- 
ian radiation spectrum everywhere. Although this is a 
fair assumption in opaque regions where the number of 
mean-free-paths is large, it fails in optically thin regions. 
Since thermalization softens the radiation spectrum, the 
assumed Planck spectrum is likely to under-predict the 
heating rate. 

Since the net affect of the approximations is somewhat 

unclear, comparison with more sophisticated radiative 
treatments would be ideal. However, there have been no 
2D or 3D non-gray simulations of low-mass star forma- 
tion. To date, the most thorough investigation of proto- 
stellar formation is presented by Masunaga & Inutsuka 
(2000). These spherically symmetric simulations follow 
the formation of 0.8 and 1.0 Mq protostars. At radii of 60 
AU, they find temperatures ranging from 20-250 K dur- 
ing the main accretion phase, while we find Tmax ~ 90 K. 
Their maximum protostellar luminosity is 25 Lq, which 
is entirely due to accretion. A few of the protostars in 
the RT simulation have higher masses and higher maxi- 
mum luminosities, but the gas temperature distributions 
on average appear similar (see Figure 4). However, the 
disparity in maximum temperature may be attributable 
to either differences in the radiation schemes or initial 
conditions and geometry. Future work will investigate 
the effects of 3D multi-frequency radiative transfer on 
low-mass star formation. 

5. CONCLUSIONS 

We perform gravito-radiation-hydrodynamic simula- 
tions to explore the effect of radiation feedback on the 
process of low-mass star formation. We compare our 
calculation with a similar one using an approximately 
isothermal equation of state in lieu of radiative trans- 
fer. We find that the inclusion of radiative transfer has 
a profound effect on the gas temperature distribution, 
accretion, and final stellar masses. 

We confirm the finding of Bate (2009b) that additional 
heating provided by radiative transfer stabilizes proto- 
stellar disks and suppresses small scale fragmentation 
that would otherwise result in brown dwarfs. However, 
we also find that the vast majority of the heating comes 
from protostellar radiation, rather than from compres- 
sion or viscous dissipation. Thus calculations that ne- 
glect radiative feedback from protostars, either because 
they use approximations for radiative effects that are 
incapable of including it (e.g.. Bate et al. 2003; Clark 
et al. 2008) or because the explicitly omit it (e.g., Bate 
2009b), significantly underestimate the gas temperature 
and thus the strength of this effect. More generally, we 
find that, due to significant variations in the tempera- 
ture with time, no scheme that does not explicitly include 
time-dependent protostellar heating is able to adequately 
follow fragmentation on scales smaller than ^0.05 pc. 

We find that due to the increased thermal support in 



the protostellar disks, accretion is smoother and less vari- 
able with radiative feedback. However, we show that for 
low-mass star formation the heating is local and limited 
to the volume within the protostellar cores. As a result, 
pre-existing sources do not inhibit turbulent fragmenta- 
tion elsewhere in the domain. 

We find that the mean accretion rate increases with 
final stellar mass so that the star formation time is only 
a weak function of mass. This is inconsistent with the 
standard Shu (1977) picture, but it is qualitatively con- 
sistent with the McKee & Tan (2003) result for the tur- 
bulent core model, where the star formation time varies 
as the final stellar mass to the 1/4 power. 

The magnitude and variability of protostellar luminos- 
ity is of significant observational interest. If accretion 
contributes a substantial portion of the total luminos- 
ity emitted by young protostars, then upper limits for 
protostellar accretion rates can be inferred directly from 
the observed luminosity. This may give clues about the 
formation timescale and the accretion process while the 
protostars are deeply embedded and cannot be directly 
imaged. In a future paper we will examine the "luminos- 
ity problem" and compare with embedded Class and 
Class I protostars. 

Our larger NRT and RT simulations are performed at 
a maximum resolution of 32 AU, so it is possible that 
a few of our cores form stars that otherwise would have 
become thermally supported or turbulently disrupted in 
a higher resolution calculation. Thus, higher resolution 
calculations would be desirable for further work. Al- 
though we find that the inclusion of radiative transfer 
has a similar impact as magnetic fields on fragmenta- 
tion and accretion, simulations examining the interplay 
of magnetic fields and radiative transfer are important. 
To asses the accuracy of our radiative transfer approxi- 
mations, further simulations with multi-frequency treat- 
ment in multi-dimensions with improved dust modeling 
are also necessary. 
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APPENDIX 

A. THE STAR PARTICLE ALGORITHM 

In this appendix we describe the details of our "star particle" algorithm we use to represent protostars. Appendix 
A describes how the star particle algorithm functions within the larger ORION code, while Appendix B describes the 
protostellar evolution code that we use to determine the luminosities of our stars. This division is useful because, 
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from the standpoint of the ORION code, a star particle is characterized by only four quantities: mass, position, 
momentum, and luminosity. The luminosity is determined by the protostellar evolution model outlined in Appendix 
B that is attached to each star particle, but the only output of this model that is visible to the remainder of the code 
is luminosity. 

In a calculation using star particles, we add a set of additional steps to every update cycle on the finest AMR level, 
so that the cycle becomes 

1. Hydrodynamic update for gas 

2. Gravity update for gas 

3. Radiation update, including stellar luminosity 

4. Star particle update 

(a) Sink particle update 

(b) Stellar model update 

Steps (1) - (3) are the ordinary parts of the update that we would perform even if no star particles were present. In 
steps (1) and (2) star particles have no direct effect (since they do not interact hydro dynamically, and we handle their 
gravitational interactions with the gas in an operator split manner in step (4a) . 

In step (3), star particles act as sources of luminosity, as indicated in equation (7). We implement this numerically 
as follows: let L„ and x„ be the luminosity and position of star particle n. Our code uses the Krumholz et al. (2007b) 
radiation-hydrodynamic algorithm, in which we split the radiation quantities into those to be handled explicitly and 
those to be handled implicitly. We therefore write the evolution equation to be solved during the radiation step as 

"q^ — fe— rad "I" fi— radj (■^-'-) 

where q = (p, pv, pe, E) is the state vector describing the gas and radiation in a cell, the explicit update vector fo-rad 
is the same as in the standard Krumholz et al. algorithm (their equation 52)^, and the implicit update is modified to 
be 
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Here — x„) is a weighting function that depends on the distance between the location of the cell center x and 
the location of the star x„. The weighting function has the property that the sum of W(x — x„) over all cells is unity, 
and that VK(x — x„) = for |x — x„| larger than some specified value. For the computations we present in this paper 
we use the same weighting function as we use for accretion (equation (13) of Krumholz et al. 2004). However, we have 
experimented with other weighting functions, including truncated Gaussians, top-hats, and delta functions, and we 
find that the choice makes very little difference because radiation injected into a small volume of the computational 
grid almost immediately relaxes to a configuration determined by diffusion. With this modification to fi_rad, our 
update procedure is the same as described in Krumholz et al. (2007b). 

Step (4a) is the ordinary sink particle method of Krumholz et al. (2004), so we only summarize it here and refer 
readers to that paper for a detailed description and the results of numerous tests. We first create new particles in any 
cell whose density exceeds the Jeans density on the maximum AMR level (i.e., where equation (12) is not satisfied.) 
Next we merge star particles whose accretion zones, defined to be 4 cells in radius, overlap. This ensures that we 
combine multiple sink particles created in adjacent cells that simultaneously exceed the Jeans density, or multiple sink 
particles created in the same cell during consecutive time steps. Then we transfer mass from the computational grid 
onto existing sink particles. Accretion happens within a radius of 4 cells around each sink particle. The amount of 
mass that a sink particle accretes is determined by fitting the flow around it to a Bondi flow, reduced to account for 
an angular momentum barrier that would prevent material from reaching the computational cell in which the sink 
particle is located. The division of mass accreted among cells inside the 4-cell accretion zone is determined by a 
weighting function. The accretion process leaves the radial velocity, angular momentum, and specific internal energy 
of the gas on the computational grid imchanged (in the frame co- moving with the sink particle), and it conserves mass, 
momentum, and energy to machine precision. Next we calculate the gravitational force between every sink particle 
and the gas in every cell using a direct 1/r^ force computation (since the number of particles is small), and modify 
the momenta of the sink particles and the momenta and energies of the gas cells appropriately. Finally we update the 
positions and velocities of each sink particle under their mutual gravitational interaction, using a simple N-body code. 
Forces are again computed via direct 1/r^ sums. 

Once the sink particle update is complete, we proceed to update the protostellar evolution model that is attached 
to each star particle. 

Note that our notation here differs slightly from that of Krumholz ct al. (2007b), in that we follow the standard astrophysics convention 
in which k is the specific opacity, while Krumholz et al. (2007b) follow the radiation-hydrodyanmic convention in which k is the total 
opacity. As a result, any opacity k that appears in the Krumholz et al. (2007b) equations is replaced by Kp here. 
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B. PROTOSTELLAR EVOLUTION MODEL 

Step (4b) of the update cycle described in Appendix A involves advancing the internal state of each star particle. 

The primary purpose of this procedure is to determine the stellar luminosity for use in step (3). Wc determine the 
luminosity using a simple one-zone protostellar evolution model introduced by Nakano ct al. (1995) and extended 
by Nakano ct al. (2000) and Tan & McKcc (2004). The model has been calibrated to match the detailed numerical 
calculations of Palla & Stabler (1991, 1992), and it agrees to ~ 10%. The numerical parameters we use for the 
calculations in this paper are based on this calibration, but we note that after we began this work Hosokawa & 
Omukai (2009) published calculations suggesting that slightly different values would improve the model's accuracy. 
We recommend that Hosokawa & Omukai's values be used in future work. 

Before diving into the details of the numerical implementation, it is helpful to give an overview of the model. The 
model essentially treats the star as a polytrope whose contraction is governed by energy conservation. The star evolves 
through a series of distinct phases, which we refer to as "pre-collapse" , "no burning" , "core deuterium burning at fixed 
Tc" , "core deuterium burning at variable Tg" , "shell deuterium burning" , and "main sequence" . The "pre-collapse" 
phase corresponds to the very low mass stage (m < 0.01 M©) when the collapsed mass is not sufficient to dissociate 
H2 and produce second collapse to stellar densities (Masunaga & Inutsuka 2000). During this phase the object is not 
yet a star. "No burning" corresponds to the phase when the object has collapsed to stellar densities, but has not 
yet reached the core temperatiirc Tc ~ 1.5 x 10^' K required to ignite deuterium, and its radiation is powered purely 
by gravitational contraction. During this phase the star is imperfectly convective. The next stage, "core deuterium 
burning at fixed Tc", begins when the star ignites deuterium. While the deuterium supply lasts, core deuterium 
burning acts as a thermostat that keeps the core temperature fixed and the star fully convective. Once the deuterium 
is exhausted, the star begins the "core deuterium burning at variable Tc" phase, during which the core temperature 
continues to rise. The star remains fully convective, and new deuterium arriving on the star is rapidly dragged down to 
the center and burned. The rising core temperature reduces the star's opacity, and eventually this shuts off convection 
in the stellar core, beginning the "shell deuterium burning" phase. At the start of this phase the star changes to a 
radiative structure and its radius swells; deuterium burning continues in a shell around the radiative core. Finally the 
star contracts enough for its core temperature to reach Tc « lO'' K, at which point it ignites hydrogen and the star 
stabilizes on the main sequence, the final evolutionary phase in our model. 

In the following sections, we give the details of our numerical implementation of this model. 

Initialization and Update Cycle 

When a star is first created, its mass is always below 0.01 M© and thus in the "pre-collapse" state. We do not 

initialize our protostellar evolution model imtil the mass exceeds 0.01 M© - prior to this point star particles are 
characterized only by a mass and have zero luminosity. On the first time step when the mass exceeds 0.01 Mq, we 
change the state to "no burning" . Thereafter each star particle is characterized by a radius r, a polytropic index n, 
and a mass of gas from which deuterium has not yet been burned, m^. We initialize these quantities to 

-2.«3(..^?^)"' (BI) 



n = 5 — 3 



^lO-5A^0yr 
1.475 + 0.07 logio 
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where At and Am are the size of the time step when the star passes 0.01 M© and the amount of mass accreted during 
it. If this produces a value of n below 1.5 or greater than 3.0, we set n = 1.5 or n = 3.0. These fitting formulae are 
purely empirical calibrations designed to match Palla & Stabler (1991, 1992). The choice of n intermediate between 
1.5 and 3.0 corresponds to imperfect convection. 

Once a star particle has been initialized and its state set to "no burning", during each time step we perform the 
following operations: 

1. Update the radius and the deuterium mass 

2. Compute the new luminosity 

3. Advance to the next evolutionary phase 
We describe each of these operations below. 

Evolution of the Radius and Deuterium Mass 

Once a star reaches the "main sequence" evolutionary phase, we simply set its radius equal to the radius of a zero- 
age main sequence star of the same mass, which we compute using the fitting formula of Tout et al. (1996) for Solar 
metallicity. Before this point we treat the star as an accreting polytrope of fixed index. For such an object, in a time 
step of size At during which the star gains a mass Am, the radius changes by an amount Ar given by a discretized 
version of equation (5.8) of Nakano et al. (2000): 

„Am / l-fk 1 dlogl3\ ^At f r , , , 

Ar = 2 1 ^ + --pr^ r - 2— _ {L,^,+Li - Ln)r B4 

m \ Ggp 2diogmJ agp \Gm^ J 
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Here flg = 3/(5 — n) is the coefBcient describing the gravitational binding energy of a polytrope, fi is the mean ratio 
of the gas pressure to the total gas plus radiation pressure in the star, fk is the fraction of the kinetic energy of the 
infalling material that is radiated away in the inner accretion disk before it reaches the stellar surface, Lint is the 
luminosity leaving the stellar interior, Li is the rate at which energy must be supplied to dissociate and ionize the 
incoming gas, and .Ld is the rate at which energy is supplied by deuterium burning. 

In this equation we adopt fk = 0.5, the standard value for an a disk. For /?, the low-mass stars we discuss in this 
paper have negligible radiation pressure and so /? = 1 to very good approximation. In general, however, we determine 
(3 and d\ogj3 / dlogm by pre-computing a table of /3 values for polytropes as a function of polytropic index n and mass 
m, and then interpolating within that table. For n = 3 interpolation is unnecessary and we instead obtain (3 by solving 
the Eddington quartic 



a \iJLmn 



/3' 



/3 4 



(B5) 



where Pc and pc are the central pressure and density of the polytrope (which are also stored in a pre-computed table 
as a function of n), and = 0.613 is the mean molecular weight for fully ionized gas of Solar composition. 
For the luminosity emanating from the stellar interior we adopt 



Lint = max (Lms, -^h) , 



(B6) 



where Lms is the luminosity of a main sequence star of mass m, which we compute using the fitting formula of Tout 
et al. (1996) for Solar metallicity, and Lh = Anr'^aT^ is the luminosity of a star on the Hayashi track, with a surface 
temperature Th = 3000 K. For the luminosity required to ionize and dissociate the incoming material we use 



Lj = 2.5 LfT) 



(Arrt/Ai) 
10-5 Mq yr- 



(B7) 



which corresponds to assuming that this process requires 16.0 eV per hydrogen nucleus. The deuterium luminosity 
depends on the evolutionary stage. In the "pre-coUapse" and "no burning" phases, Ld = 0. In the "core burning at 
fixed Tc" phase, we set the deuterium luminosity to the value required to keep the central temperature at a constant 
value Tc = 1.5 x 10*^ K. This is (equation (5.13) of Nakano et al. 2000) 
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where (3c = PcksTc/ {jJiriYiPc) is the ratio of gas pressure to total pressure at the center of the polytrope. In all 
subsequent phases, deuterium is burned as quickly as it is accreted, so we take 



Ld = 15 Lq 



{Am/ At) 
10-5 AfQ yr- 



(B9) 



which corresponds to assuming an energy release of 100 eV per gram of gas, appropriate for deuterium burning in a 
gas where the deuterium abundance is D/H = 2.5 x 10"^. Finally, we update the mass of material that still contains 
deuterium simply based on Ld- The change in unburned mass is 



Amd = Am - 10~^Mg 



L 



D 



15 Lr. 



At 
yr 



(BIO) 



Com,puting the Luminosity 

From the standpoint of the rest of the code, the only quantity of any consequence is the luminosity, since this is 
what enters as a source term in step (3). The luminosity radiated away from the star consists of three parts: 



— -^int + -^disk, 



(Bll) 



where Lint is the luminosity leaving the stellar interior as defined above, Lace is the luminosity radiated outward at 
the accretion shock, and L^isk is the luminosity released by material in traversing the inner disk. These in turn are 
given by 



-^acc — /acc/fe 



GmAm/At 



-t'disk = (1 - fk) 
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(B12) 
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where fk = 0.5 as defined above, and /acc is the fraction of the accretion power that is radiated away as light rather 
than being used to drive a wind. Although we do not explicitly include a protostellar outflow in this calculation, 
we take /acc = 0.5 so that we do not overestimate the accretion luminosity by assuming that the all the accretion 
power comes out radiatively rather than mechanically. Thus, we assume a total radiative efficiency of 75%. Although 
this value is consistent with x-wind models (Ostriker & Shu 1995), neither x-wind or disk- wind models definitively 
constrain the total conversion of accretion energy into radiation, and we treat this as a free parameter. 
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Advancing the Evolutionary State 

The final pieces of our protostellar evolution model are the rules for determining when to change the evolutionary 
state, and for determining what happens at such a change. Our rules are as follows: if the current state is "no burning" , 
then at the end of each time step we compute the central temperature by numerically solving the equation 

Pc = + laT^, (B14) 

limn 3 

where Pc and pc are determined from the current mass, radius, and polytropic index. If Tc > 1.5 x 10^ K, we change 
the evolutionary state to "core burning at fixed Tc" and we change the polytropic index to n = 1.5. 

If the current evolutionary state is "core burning at fixed Tc" , then we check to make sure that uid > after we 
update the unburned deuterium mass with equation (BIO). If not, then the deuterium has been exhausted and we 
change the state to "core burning at variable Tc" . 

If the current state is "core burning at variable Tc" , we decide whether a radiative zone has formed by comparing 
the luminosity being generated by deuterium burning, Ld, to the hmiinosity of a zero-age main sequence star of the 
same mass, Lms. We switch the state to "shell deuterium burning" when Lo/L^is > ./rad = 0.33. At this point we 
also change the polytropic index to n = 3 and increase the radius by a factor of 2.1, representing a swelling of the star 
due to formation of the radiative barrier. 

Finally, if the state is "shell burning", we compare the radius r at the end of every time step to the radius of a 
zero-age main sequence star of the same mass. Once the radius reaches the main sequence radius, we switch the state 
to "main sequence" , our final evolutionary state. 
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